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Abstract 

A soft ellipsoid contact potential model for a pair of biaxial ellipsoidal molecules is proposed 
which considers the configuration dependent energy anisotropy explicitly along with their geo¬ 
metrical aspects. We performed Molecular Dynamics simulation study to generate both biaxial 
smectic and nematic phases using this new potential. 
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Liquid crystal forming compounds are complex organic molecules, having in general, cylindri- 
cally symmetric rod-like or disc-like structures,. Theoretical models, developed to study phase 
transitions, structures and dynamics of these systems, essentially strive to relate important molec¬ 
ular interactions to the experimentally observed phase behaviour and properties. In this respect, 
computer simulation is one of the most efficient methods which considers detailed microscopic in¬ 
teractions theoretically to explain and analyze macroscopic experimental observations. Non-polar 
liquid crystal molecules interact both through short-range hard interactions and long range attrac¬ 
tive interactions originated primarily from electrostatic interactions [I]. To study a liquid crystal 
system using computer simulation two major steps are to be followed. The first one is the modelling 
of a suitable potential which takes into account the basic electrostatic interactions responsible for 
the generation of one (or more) specific phase(s) and the second step is to use an appropriate sim¬ 
ulation method. These models may be broadly classified as full atomic, united atom, site-site and 
single site interaction models. The first three class of models are not computationally efficient as 
there are multi-site interaction terms present in the potential. The fourth type of model potentials 
can be used efficiently for computer simulation studies of these complex molecular systems in a 
realistic way. 

Single-site potentials require determination of a, distance of closest approach as a function 
of orientation and center of mass distance for a pair of molecules and it is a difficult task to 
estimate this correctly in case of anisotropic molecular systems like liquid crystals. Additionally, 
the potentials must capture both the geometric and energy anisotropies of the systems consistently. 
To develop single-site model potentials which are accurate, relatively simple and computationally 
efficient, a rod-like rigid molecule may be approximated as a prolate ellipsoid. The main idea behind 
it is the representation of the interactions between a pair of liquid crystal molecules by their joint 
probability, where each molecule has Gaussian charge distribution centred around the molecular 
centroids. Corner [2]developed a model potential based on this idea which makes the analytical 
calculations remarkably simplified. Further significant development of this class of potential was 
done by Berne’s group which led to widely used Gay-Berne(G-B) potential [3]. These models 
were basically modified form of Lennard-Jones potential where the energy strength parameter e 
and energy range parameter a depended on molecular orientations and positions. In this form 
it provided both attractive and repulsive contributions of the intermolecular interactions of rigid, 
uniaxial molecules, in a computationally tractable way. These potentials known as Gaussian overlap 
potentials (GOP) have become popular in simulation studies due to their comparatively simplihed 
approach in determining the energy contribution of a pair of ellipsoidal molecules. However, the G-B 
potential in its original form being generally restricted to uniaxial molecules is not an effective model 
when the multiplicity of phases adopted by liquid crystal molecules are to be studied. Modifications 
[1] and some parametrizations [5] of Gay-Berne potential have been proposed in order to study the 
behaviour of biaxial phases but these GOP models are not in conformity with the actual geometry 
of the molecule. The range function and the energy strength function in these potentials do not 
correctly estimate the interactions between the non-uniaxial ellipsoidal particles. Most of the real 
liquid crystal molecules do not have cylindrical uniaxial structures, therefore, biaxial contribution 
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to the potential is essential when one studies various phase structures adopted by the real molecules. 
Though much progress in the experimental study have been achieved, there has been comparatively 
little theoretical work done on this important class of liquid crystal forming molecules which have 
biaxial structures. 

Perram-Wertheim [6] ellipsoid contact potential (ECP) had advantages over Gaussian overlap 
type potential because the ECP calculated geometrical contact point accurately for ellipsoidal 
structures including biaxial molecules. ECP represented the ellipsoids as quadratic forms to solve 
the contact function accurately by doing an optimization calculation. Despite this fact, ECP has 
not yet been widely used to study liquid crystals because it did not differentiate in energies due to 
different molecular configurations, thereby neglecting the proper weightage of the energy strength 
contribution to the potential. Energy strength parameter e was taken as a constant in the ECP 
therefore well depth in this attractive-repulsive potential did not show correct behaviour. In our 
earlier work [8], we modelled a soft ellipsoid contact potential for uniaxial ellipsoidal molecules, 
considering the configuration dependent energy anisotropy along with their geometrical aspects. In 
this paper we generalize the earlier model by including biaxial contribution in the pair interaction 
potential for ellipsoidal molecules. Explicit analytical expressions for proper range and strength 
parameters alongwith related forces and torques are also provided. In addition, the molecular 
dynamics simulation study shows that this new soft ellipsoidal pair potential for biaxial molecules 
is suitable for generating both biaxial and uniaxial liquid crystal phases successfully in a single 
component rod-like molecular system. 

In our model the charge distribution of molecules are described as ellipsoidal and are approx¬ 
imately represented by 3-D gaussian functions, i.e. the density of the centre of force along the 
axis of a molecule is a Gaussian function of distance from the center of the ellipsoid. The total 
interaction of two molecules is again a Gaussian distribution and is a function of their intercenter 
distance and relative orientations. 

The pair interaction potential of two ellipsoidal molecules can be evaluated from the joint proba¬ 
bility distribution of two Gaussians following convolution integral principle. This joint distribution 
function is proportional to 

exp[r'''.H~^.r] 

where r is the intermolecular vector connecting the centers of the ellipsoids and the matrix H 
depends on the axes lengths and the relative orientations of three principal semi-axes of the two 
interacting ellipsoids. Both the energy range function a i.e. the permitted closest distance between 
molecular pairs and the energy strength function e which gives estimation of potential depths and 
widths, depend on this form of joint probability distribution functions. 

The model pair potential which is a modified version of Lennard-Jones type interaction potential 
for the prolate ellipsoidal molecules can be represented in a shifted form as 
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( 1 ) 


U = 4e(Q;,/3,f)[(- 

r — a 

_oo_ 

r - a{a,f3,r) + ao 

where, a and (3, are the orthogonal rotation matrices for transformation from space-fixed to 
molecule fixed coordinates for molecules 1 and 2 respectively and the value of ao is equal to the 
smallest semi-axis length of the ellipsoid. 

Unlike constant e parameter in ECP model, our model considers relative configuration depen¬ 
dent e for a molecule pair because they are characterzised by highly non-spherical charge distribu¬ 
tions. The energy strength function e in our model potential is taken as 

e = eoei(d,/3)e^(Q;,/3,f). (2) 

Here, cq is a constant term and fi, v are adjusting parameters for relative well depth variation in 
different compounds; ei is a function of a and /3; 62 depends additionally on f, which is the unit 
vector along the intermolecular vector connecting the centers of the ellipsoids. 


'^0 _^12 

{a,j3,f) + fJo 


To calculate 62 , we consider diagonal energy matrix for a biaxial molecule as 


/ 



7 = 


V 


0 
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(3) 


where €x,€y, are equal to the well depths for side by side (s-s), width to width (w-w) and end to 
end (e-e) configurations and for the biaxial model we take ex > Sy > Sz- We also take eo = ex- 
Defining matrices A and B which correspond to the first and second molecule respectively as 


A = d; B = 'y /3. 


(4) 


We write 

He = AeA -|- (1 — AejB. (5) 

where Ae is an energy related adjustable parameter such that Ag G [ 0 , 1 ] and this makes the object 
function 

Se{^e) = Ae(1 — ^e)^^ ■ He“^ • r. (6) 

an optimum. This object function gives correct estimation of the energy strength function by 
properly incorporating asymmetric configurations of a biaxial molecule pair. 

The f dependent part of the energy strength function €2 in our model potential is taken as 


e2(d, f) = Ae(1 - Xe)^^ ■ [(AeA (1 - Ae)B] ^ • f 

= Ae (1 — Ae)^''' • He“^ • r 


( 7 ) 
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We take ei for biaxial molecules as [7j 


ei(a,/3) = Ui|HEr^/^ (8) 

where = ^(A + B) and at = Y^2^fTQ[(^) + ax, ay and az are the values of the range 

function for s-s, w-w and e-e configurations respectively. We calculate a following ECP [6]. 
Defining matrices M and N related to the hrst and the second molecule respectively as 

M = a; N = /3 (9) 

where diagonal shape matrix S takes the folowing form 


/ 


0 

0 ^ 



0 

ay 

0 

( 10 ) 

V 

0 

0 

fT. ) 



and Hd = AdM + (1 — A£))N. We express the distance of closest approach for biaxial molecules a 
in the following form 

a~'^ = Xd{1 - (11) 

where Ad is the shape related adjustable parameter. For our biaxial model ax < ay < az- The 
potential is not isotropic at large distance which gives small thermodynamic contribution and 
therefore neglected for simplicity. The forces and torques are derived analytically for the biaxial 
ellipsoid. The expressions for force / and torque r are 

/ = [(3r^(2p~^^ - p~'^)r + 3o-^Ad(1 - Ad) 

aor^ 

{2p~^^ - /9“’^){k - (f.k)f} - Ae(1 - Ad) 

^62 - P”®){kD - {r.kE)f}] 

( 12 ) 


Tist mol. = ^^^^3^[-3(T^Ad(1 - AD)(k.M X k) 

{2p~^^ — p~'^)jao + Ad( 1 — AD)(kD.A x kD) 
{p-^^ - p-^)] + i^Ue:B-^M 


(e in equation (12) is the Levi-Civita tensor) 

To check the suitability of the proposed model, we performed a sample NVE Molecular Dynamics 
(MD) simulation considering a system of prolate biaxial ellipsoids, interacting through the model 
potential described in the previous section. In this simulation we take two systems having different 
aspect ratios and well depth ratios. The systems consisting of 1372 ellipsoids interacting through 
this new potential all have /r = 1, d = 2. The other set of parameters are taken as follows: 
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The ratios ax ■ CTy : Uz have values 1 ; 1.5 ; 4.5 and 1:2: 4.5 and corresponding 1 * 3 ,, I*y^ have the 
set of values (1.125,1.062, 0.16); (1.2125,1.062, 0.25) for two systems respectively. All the systems 
have exl^z = 1/30 but having different e^/e^ values 1/7,1/8.5 for two systems respectively. 

In the NVE MD simulation, 1372 molecules were confined in the cubic box obeying minimum 
image convention. The run started from a density p* = paQ = 0.01 with ellipsoidal molecules loacted 
on the sites of FCC lattice and having parallel orientation. For MD initial step translational and 
angular velocities were assigned from the Gaussian distribution of velocities at a higher temperature. 
This structure was not a stable structure at this density and it was melted at a reduced temperature 
T* = ksT/eo = 4.0 . We used this isotropic configuration which was both orientationally and 
translationally disordered, as the initial configuration. Forces and torques were calculated from 
equations (12) and (13) using leap-frog algorithm [9] . To get the most ordered structure for 
the system, we increased the density from p* = 0.01 to the desired values of p* = 0.19 and 0.16 
respectively for two set of systems with an increament size of 0.001 upto p* = 0.1 and 0.01 for the 
rest. Temperature was then lowered in steps to help the molecules set with more order. The systems 
were cooled down to a very low reduced temperature T* = ksT/eo = 0.4. The interaction potential 
was taken equal to half the simulation box length. The time step size 6t* = 6t/{ma^/eo)^^'^ = 0.0012 
during cooling process. The orientations of molecules were described by quaternions instead of 
Fulerian angles to get the singularity-free orientational equations of motion. 

To investigate the phase behaviour we started increasing the temperature T* gradually. For 
each temperature we calculated the orientational order parameters, correlation functions to identify 
the particular liquid crystal phase. To achieve an equilibrated higher temperature configuration 
from more ordered low temperature phase 2 x 10 ^ cycles were required far from transition and 
5 X 10^ near transition. 

The orientational order parameter for uniaxial phase was calculated from the largest eigen value 
obtained by diagonalization of the order parameter tensor 

Qap = - 6ap) a,l3 = x,y,z. (14) 

where was the a th component of the unit vector e* along the symmetry axis of the i 
th molecule. Corresponding eigenvector gave the director which defines the average direction of 
molecular alignment. 

The biaxial order parameter is (ii^ 2 ) = /^) ^o® ^7 ~ cos (3 sin 2 q; sin 27 ). 

Additionally, to detect phase behaviour more accurately, some structural quantities were cal¬ 
culated. Those were radial distribution function g{r) =< 6{r — rij) >ij /Airr'^p; density along the 
director g{z) =< 5{z — Zij) >ij /irR'^p, where Zij = rijcosf^ri^ was measured in the director frame 
and R is the radius of the cylindrical sampling region. 

The snapshots of the simulation study show various phases formed at different temperatures. 
At T* = 0.4 the smectic phase shows biaxiality in both cases as seen in the structures for the first 
system in Fig.l and Fig.2 and for the second system in Fig.3 and Fig.4. At T* = 1.0 it shows biaxial 
nematic phases in Fig.5 and Fig .6 for the first system and fig.7 and Fig .8 for the second system. 
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At T* > 5.0 the systems becomes isotropic. The radial distribution functions g{r) (Fig.9, Fig.ll) 
and density projection to the director g{z) (Fig. 10, Fig. 12) are plotted at different temperatures 
for both systems. The values of orientational order parameters r/ and R 22 are plotted in Fig. 13 
and Fig. 14. Our results show that the new model potential can successfully reproduce both biaxial 
smectic and biaxial nematic phases. To demonstrate this new model’s ability to generate biaxial 
liquid crystal phases we are reporting here some preliminary results. We plan to communicate 
details of phase transitional behaviour of this system in future work. 

Liquid crystals show a rich multiplicity of phases which comes from their complex molecular 
structures giving rise to complex electromagnetic interactions. Other molecules of biological impor¬ 
tance e.g. base pairs of DNA, lipids etc. basically have anisotropic biaxial structures. It is a difficult 
task to obtain the distance of closest approach as a function of orientation for these molecules. Ad¬ 
ditionally, this potential must capture the energetic strength anisotropy in a consistent way. GOP 
is a widely used model potential for the description of anisotropic behaviour of liquid crystals, 
however, it does not give a proper geometric interpretation. It can calculate the proper range and 
strength parameters for some specific configurations of identical uniaxial molecules. Perram and 
Wertheim took a different approach by introducing the ellipsoid contact function which took care of 
its true geometrical aspect by doing correct estimation of position of contact. Although applicable 
to any mixture in geometrical context it did not distinguish in energy functions between different 
relative orientations and positions of the non-spherical molecular pair. As this energy function is 
the basis of well depth estimation, it provides necessary and sufficient weightage while determining 
proper magnitude of the attaractive part of the pair potential and thus acts as the most impor¬ 
tant contributor in fine tuninng of the potential so that it generates a specific phase structure. In 
this work, considering molecular orientations and positions explicitly in both structural and energy 
strength related terms, we modelled a generalized pair potential for biaxial molecules. We further 
studied liquid crystal phases by doing molecular dymanics simulation utilizing this model potential. 
Though biaxial phases have been synthesized and analysed recently, theoretical studies including 
simulations have been done for more than three decades to understand the molecular properties 
needed for stabilization of these amazing phases. However, underlying intermolecular interactions 
responsible for the generation of these phase are still not understood properly and the present 
model potential will be able to provide essential informations about the characteristic of single 
component rod-like molecular system crucial for obtaining biaxial phases. The simulation results 
show that this new model potential reported in this paper is qualitatively suitable to reproduce 
complex mesogenic behaviour of organic biaxial molecules in a realistic way. This new potential 
has been generalized to study systems comprising of multicomponent dissimilar molecules, results 
of which will be communicated in future publication. 

The graphics software QMGA [TO] was used for producing snapshots of the phases. This work 
was supported by UGG-UPE scheme of the University of Galcutta. 
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Figure 1: Snapshot at T* = 0.4 for Cx : <7^ : = 1 : 1.5 : 4.5 



Figure 2: Snapshot (rotated) at T* = 0.4 for (Tx : <7 ^ : cr^ = 1 : 1.5 : 4.5 
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Figure 3: Snapshot at T* = 0.4 for ax '■ cry '■ cTz = ^ 2-0 : 4.5 



Figure 4: Snapshot (rotated) at T* = 0.4 for Ua; : Uy : = 1 : 2.0 : 4.5 
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Figure 5: Snapshot at T* = 1.0 for Ua: : cr^ : Uz = 1 : 1.5 : 4.5 



Figure 6: Snapshot (rotated) at T* = 1.0 for Ua; : Uy : Cz = 1 : 1.5 : 4.5 
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Figure 7: Snapshot (rotated) at T* = 1.0 for ax '■ CFy : cjz = 1 '■ 2.0 : 4.5 



Figure 8: Snapshot (rotated) at T* = 1.0 for cja; : Uy : = 1 : 2.0 : 4.5 
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Figure 9: Radial distribution function g(r) for (Tx '■ (Ty : = 1 '■ 1-5 : 4.5 



0 0.2 0.4 


z 

Figure 10: Denity projection with respect to the director g(z) for '■ (Ty : az = 1 : 1.5 : 4.5 









Figure 11: Radial distribution function g(r) for ax ■ (Jy '■ a^ 
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Figure 12: Denity projection with respect to the director g(z) for ax ■ ay : a. 














